Hands-on_Ex02

Installing and Loading the R packages

pacman::p_load(sf, raster, spatstat, tmap, tidyverse)

Spatial Data Wrangling

Importing the spatial data

childcare_sf <- st_read("/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial/child-care-services-geojson.geojson") %>%
  st_transform(crs = 3414)
Reading layer `child-care-services-geojson' from data source 
  `/Users/sharon/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial/child-care-services-geojson.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
sg_sf <- st_read(dsn = "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial/", layer="CostalOutline")
Reading layer `CostalOutline' from data source 
  `/Users/sharon/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 67 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6091 ymin: 1.16639 xmax: 104.0858 ymax: 1.471388
Geodetic CRS:  WGS 84
mpsz_sf <- st_read(dsn = "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial/", 
                layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `/Users/sharon/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

Mapping the geospatial data sets

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare_sf)+
  tm_dots()
tmap_mode('plot')
tmap mode set to plotting

Geospatial Data wrangling

childcare <- as_Spatial(childcare_sf)
mpsz <- as_Spatial(mpsz_sf)
sg <- as_Spatial(sg_sf)
childcare
class       : SpatialPointsDataFrame 
features    : 1925 
extent      : 11810.03, 45404.24, 25596.33, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 2
names       :    Name,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    Description 
min values  :   kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>100044</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>44, TELOK BLANGAH DRIVE, #01 - 19/51, SINGAPORE 100044</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td></td> </tr><tr bgcolor=""> <th>NAME</th> <td>PCF SPARKLETOTS PRESCHOOL @ TELOK BLANGAH BLK 44 (CC)</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>349C54F201805938</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20211201093837</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
max values  : kml_999,                                            <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>99982</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>35, ALLANBROOKE ROAD, SINGAPORE 099982</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td></td> </tr><tr bgcolor=""> <th>NAME</th> <td>ISLANDER PRE-SCHOOL PTE LTD</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>4F63ACF93EFABE7F</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20211201093837</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
mpsz
class       : SpatialPolygonsDataFrame 
features    : 323 
extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs 
variables   : 15
names       : OBJECTID, SUBZONE_NO, SUBZONE_N, SUBZONE_C, CA_IND, PLN_AREA_N, PLN_AREA_C,       REGION_N, REGION_C,          INC_CRC, FMEL_UPD_D,     X_ADDR,     Y_ADDR,    SHAPE_Leng,    SHAPE_Area 
min values  :        1,          1, ADMIRALTY,    AMSZ01,      N, ANG MO KIO,         AM, CENTRAL REGION,       CR, 00F5E30B5C9B7AD8,      16409,  5092.8949,  19579.069, 871.554887798, 39437.9352703 
max values  :      323,         17,    YUNNAN,    YSSZ09,      Y,     YISHUN,         YS,    WEST REGION,       WR, FFCCF172717C2EAF,      16409, 50424.7923, 49552.7904, 68083.9364708,  69748298.792 
sg
class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 103.6091, 104.0858, 1.16639, 1.471388  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 67
names       : ID_0, ISO, NAME_ENGLI,  NAME_ISO,  NAME_FAO, NAME_LOCAL, NAME_OBSOL, NAME_VARIA, NAME_NONLA, NAME_FRENC, NAME_SPANI,                 NAME_RUSSI,                      NAME_ARABI,               NAME_CHINE, WASPARTOF, ... 
value       :  205, SGP,  Singapore, SINGAPORE, Singapore,  Singapore,         NA,         NA,         NA,  Singapour,   Singapur, Сингапур, سنغافورة, 新加坡,        NA, ... 

Converting the Spatial* class into generic sp format

childcare_sp <- as(childcare, "SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")
childcare_sp
class       : SpatialPoints 
features    : 1925 
extent      : 11810.03, 45404.24, 25596.33, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
sg_sp
class       : SpatialPolygons 
features    : 1 
extent      : 103.6091, 104.0858, 1.16639, 1.471388  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 

Converting the generic sp format into spatstat’s ppp format

childcare_ppp <- as.ppp(childcare_sf)
Warning in as.ppp.sf(childcare_sf): only first attribute column is used for
marks
childcare_ppp
Marked planar point pattern: 1925 points
marks are of storage type  'character'
window: rectangle = [11810.03, 45404.24] x [25596.33, 49300.88] units
plot(childcare_ppp)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 1925 symbols are shown in the symbol map

summary(childcare_ppp)
Marked planar point pattern:  1925 points
Average intensity 2.417323e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1925 character character 

Window: rectangle = [11810.03, 45404.24] x [25596.33, 49300.88] units
                    (33590 x 23700 units)
Window area = 796335000 square units

Handling duplicated points

any(duplicated(childcare_ppp))
[1] FALSE
multiplicity(childcare_ppp)
   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [445] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [482] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [519] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [556] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [593] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [630] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [667] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [704] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [741] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [778] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [815] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [852] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [889] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [926] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [963] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1000] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1037] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1074] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1111] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1148] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1185] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1222] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1259] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1296] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1333] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1370] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1407] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1444] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1481] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1518] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1555] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1592] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1629] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1666] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1703] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1740] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1777] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1814] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1851] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1888] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1925] 1
sum(multiplicity(childcare_ppp) > 1)
[1] 0
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare) +
  tm_dots(alpha=0.4, 
          size=0.05)
tmap_mode('plot')
tmap mode set to plotting
childcare_ppp_jit <- rjitter(childcare_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)
any(duplicated(childcare_ppp_jit))
[1] FALSE
# Transform the CRS to SVY21 (EPSG:3414)
sg_sf_projected <- st_transform(sg_sf, 3414)

# Check if the CRS transformation was successful
st_crs(sg_sf_projected)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

Creating owin object

sg_owin <- as.owin(sg_sf_projected)
plot(sg_owin)

summary(sg_owin)
Window: polygonal boundary
34 separate polygons (no holes)
            vertices        area relative.area
polygon 1       4033 604202000.0      8.67e-01
polygon 2         38    414592.0      5.95e-04
polygon 3         16     19671.9      2.82e-05
polygon 4        537  24890800.0      3.57e-02
polygon 5        106   1949160.0      2.80e-03
polygon 6        261  10256000.0      1.47e-02
polygon 7         36     38729.5      5.56e-05
polygon 8         76   1087280.0      1.56e-03
polygon 9         62    920917.0      1.32e-03
polygon 10        96    573783.0      8.23e-04
polygon 11        34    126330.0      1.81e-04
polygon 12        13     20190.5      2.90e-05
polygon 13        18     34466.7      4.94e-05
polygon 14        37    239043.0      3.43e-04
polygon 15       320   4789260.0      6.87e-03
polygon 16       473  26400400.0      3.79e-02
polygon 17        21     47547.0      6.82e-05
polygon 18        77   1430070.0      2.05e-03
polygon 19       192   4796520.0      6.88e-03
polygon 20       132   2985250.0      4.28e-03
polygon 21        75    955310.0      1.37e-03
polygon 22        19     63805.9      9.15e-05
polygon 23       105    832107.0      1.19e-03
polygon 24        21     35903.4      5.15e-05
polygon 25        43    139236.0      2.00e-04
polygon 26        23     71722.4      1.03e-04
polygon 27        70    483000.0      6.93e-04
polygon 28       207   3885250.0      5.57e-03
polygon 29        16     35110.5      5.04e-05
polygon 30        26     72188.5      1.04e-04
polygon 31       156   2364690.0      3.39e-03
polygon 32        67    732165.0      1.05e-03
polygon 33       104   1100540.0      1.58e-03
polygon 34        82   1034550.0      1.48e-03
enclosing rectangle: [3040.59, 56097.76] x [16599.19, 50324.13] units
                     (53060 x 33720 units)
Window area = 697027000 square units
Fraction of frame area: 0.39

Combining point events object and owin object

childcareSG_ppp = childcare_ppp[sg_owin]
summary(childcareSG_ppp)
Marked planar point pattern:  1924 points
Average intensity 2.760294e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1924 character character 

Window: polygonal boundary
34 separate polygons (no holes)
            vertices        area relative.area
polygon 1       4033 604202000.0      8.67e-01
polygon 2         38    414592.0      5.95e-04
polygon 3         16     19671.9      2.82e-05
polygon 4        537  24890800.0      3.57e-02
polygon 5        106   1949160.0      2.80e-03
polygon 6        261  10256000.0      1.47e-02
polygon 7         36     38729.5      5.56e-05
polygon 8         76   1087280.0      1.56e-03
polygon 9         62    920917.0      1.32e-03
polygon 10        96    573783.0      8.23e-04
polygon 11        34    126330.0      1.81e-04
polygon 12        13     20190.5      2.90e-05
polygon 13        18     34466.7      4.94e-05
polygon 14        37    239043.0      3.43e-04
polygon 15       320   4789260.0      6.87e-03
polygon 16       473  26400400.0      3.79e-02
polygon 17        21     47547.0      6.82e-05
polygon 18        77   1430070.0      2.05e-03
polygon 19       192   4796520.0      6.88e-03
polygon 20       132   2985250.0      4.28e-03
polygon 21        75    955310.0      1.37e-03
polygon 22        19     63805.9      9.15e-05
polygon 23       105    832107.0      1.19e-03
polygon 24        21     35903.4      5.15e-05
polygon 25        43    139236.0      2.00e-04
polygon 26        23     71722.4      1.03e-04
polygon 27        70    483000.0      6.93e-04
polygon 28       207   3885250.0      5.57e-03
polygon 29        16     35110.5      5.04e-05
polygon 30        26     72188.5      1.04e-04
polygon 31       156   2364690.0      3.39e-03
polygon 32        67    732165.0      1.05e-03
polygon 33       104   1100540.0      1.58e-03
polygon 34        82   1034550.0      1.48e-03
enclosing rectangle: [3040.59, 56097.76] x [16599.19, 50324.13] units
                     (53060 x 33720 units)
Window area = 697027000 square units
Fraction of frame area: 0.39

First-order Spatial Point Patterns Analysis

kde_childcareSG_bw <- density(childcareSG_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 
plot(kde_childcareSG_bw)

bw <- bw.diggle(childcareSG_ppp)
bw
   sigma 
305.2404 
childcareSG_ppp.km <- rescale.ppp(childcareSG_ppp, 1000, "km")
childcareSG_ppp.km <- rescale.ppp(childcareSG_ppp, 1000, "km")
kde_childcareSG.bw <- density(childcareSG_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG.bw)

Working with different automatic badwidth methods

bw.CvL(childcareSG_ppp.km)
  sigma 
4.53932 
bw.scott(childcareSG_ppp.km)
 sigma.x  sigma.y 
2.160434 1.395302 
bw.ppl(childcareSG_ppp.km)
   sigma 
0.389485 
bw.diggle(childcareSG_ppp.km)
    sigma 
0.3052404 
kde_childcareSG.ppl <- density(childcareSG_ppp.km, 
                               sigma=bw.ppl, 
                               edge=TRUE,
                               kernel="gaussian")
par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "bw.diggle")
plot(kde_childcareSG.ppl, main = "bw.ppl")

Working with different kernel methods

par(mfrow=c(2,2))
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="gaussian"), 
     main="Gaussian")
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="epanechnikov"), 
     main="Epanechnikov")
Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="quartic"), 
     main="Quartic")
Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="disc"), 
     main="Disc")
Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel

Fixed and Adaptive KDE

kde_childcareSG_600 <- density(childcareSG_ppp.km, sigma=0.6, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG_600)

kde_childcareSG_adaptive <- adaptive.density(childcareSG_ppp.km, method="kernel")
plot(kde_childcareSG_adaptive)

par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "Fixed bandwidth")
plot(kde_childcareSG_adaptive, main = "Adaptive bandwidth")

# Load the required libraries
library(spatstat.geom)  
library(sp)         

# Step 1: Extract pixel values (z) and grid coordinates (x, y)
x_coords <- kde_childcareSG.bw$xcol   # x coordinates of the grid
y_coords <- kde_childcareSG.bw$yrow   # y coordinates of the grid
z_values <- t(kde_childcareSG.bw$v)   # Pixel values (transposed to match coordinate orientation)

# Step 2: Create a grid of the coordinates
grid <- expand.grid(x = x_coords, y = y_coords)

# Step 3: Create a SpatialPixelsDataFrame (used for gridded data)
sp_grid <- SpatialPixelsDataFrame(points = grid, 
                                  data = data.frame(value = as.vector(z_values)),
                                  proj4string = CRS("+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs"))

# Step 4: Convert to SpatialGridDataFrame
sp_grid_df <- as(sp_grid, "SpatialGridDataFrame")

# Check the structure
summary(sp_grid_df)
Object of class SpatialGridDataFrame
Coordinates:
        min      max
x  3.040593 56.09776
y 16.599186 50.32413
Is projected: TRUE 
proj4string :
[+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1
+x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs]
Grid attributes:
  cellcentre.offset  cellsize cells.dim
x          3.247848 0.4145091       128
y         16.730924 0.2634761       128
Data attributes:
     value       
 Min.   : 0.000  
 1st Qu.: 0.000  
 Median : 0.109  
 Mean   : 2.767  
 3rd Qu.: 4.598  
 Max.   :35.129  
 NA's   :9998    
#gridded_kde_childcareSG_bw <- as.SpatialGridDataFrame.im(kde_childcareSG.bw)
#spplot(gridded_kde_childcareSG_bw)
# Load the raster package
library(raster)
kde_childcareSG_bw_raster <- raster(kde_childcareSG.bw)
kde_childcareSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4145091, 0.2634761  (x, y)
extent     : 3.040593, 56.09776, 16.59919, 50.32413  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : -7.060279e-15, 35.12852  (min, max)
projection(kde_childcareSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4145091, 0.2634761  (x, y)
extent     : 3.040593, 56.09776, 16.59919, 50.32413  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : layer 
values     : -7.060279e-15, 35.12852  (min, max)

Visualising the output in tmap

# Load the tmap package
library(tmap)

tmap_mode("plot")  # For interactive maps
tmap mode set to plotting
tm_shape(kde_childcareSG_bw_raster) + 
  tm_raster("layer", palette = "viridis") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE)

Comparing Spatial Point Patterns using KDE

# Load the dplyr package
library(dplyr)
pg <- mpsz_sf %>%
  filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_sf %>%
  filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_sf %>%
  filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_sf %>%
  filter(PLN_AREA_N == "JURONG WEST")
pg <- mpsz_sf %>%
  filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_sf %>%
  filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_sf %>%
  filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_sf %>%
  filter(PLN_AREA_N == "JURONG WEST")
par(mfrow=c(2,2))
plot(pg, main = "Ponggol")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

plot(tm, main = "Tampines")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

plot(ck, main = "Choa Chu Kang")
Warning: plotting the first 10 out of 15 attributes; use max.plot = 15 to plot
all

plot(jw, main = "Jurong West")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

# Load the spatstat.geom package
library(spatstat.geom)
pg_owin = as.owin(pg)
tm_owin = as.owin(tm)
ck_owin = as.owin(ck)
jw_owin = as.owin(jw)
childcare_pg_ppp = childcare_ppp_jit[pg_owin]
childcare_tm_ppp = childcare_ppp_jit[tm_owin]
childcare_ck_ppp = childcare_ppp_jit[ck_owin]
childcare_jw_ppp = childcare_ppp_jit[jw_owin]
# Load the spatstat package
library(spatstat)
childcare_pg_ppp.km <- rescale.ppp(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp.km <- rescale.ppp(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp.km <- rescale.ppp(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp.km <- rescale.ppp(childcare_jw_ppp, 1000, "km")

# Calculate bandwidths using bw.ppl()
bw_pg <- bw.ppl(childcare_pg_ppp.km)
bw_tm <- bw.ppl(childcare_tm_ppp.km)
bw_ck <- bw.ppl(childcare_ck_ppp.km)
bw_jw <- bw.ppl(childcare_jw_ppp.km)

# Perform kernel density estimation with the calculated bandwidth
kde_pg <- density(childcare_pg_ppp.km, sigma = bw_pg, edge = TRUE, kernel = "gaussian")
kde_tm <- density(childcare_tm_ppp.km, sigma = bw_tm, edge = TRUE, kernel = "gaussian")
kde_ck <- density(childcare_ck_ppp.km, sigma = bw_ck, edge = TRUE, kernel = "gaussian")
kde_jw <- density(childcare_jw_ppp.km, sigma = bw_jw, edge = TRUE, kernel = "gaussian")


# Plot the kernel density estimates
par(mfrow = c(2, 2))
plot(kde_pg, main = "Punggol")
plot(kde_tm, main = "Tampines")
plot(kde_ck, main = "Choa Chu Kang")
plot(kde_jw, main = "Jurong West")

par(mfrow=c(2,2))
plot(density(childcare_pg_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")
plot(density(childcare_tm_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tempines")
plot(density(childcare_ck_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Choa Chu Kang")
plot(density(childcare_jw_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="JUrong West")

par(mfrow=c(2,2))
plot(density(childcare_ck_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Chou Chu Kang")
plot(density(childcare_jw_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="JUrong West")
plot(density(childcare_pg_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")
plot(density(childcare_tm_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tampines")

Nearest Neighbour Analysis

clarkevans.test(childcareSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("clustered"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcareSG_ppp
R = 0.52559, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
clarkevans.test(childcare_ck_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_ck_ppp
R = 0.88703, p-value = 0.063
alternative hypothesis: two-sided
clarkevans.test(childcare_tm_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_tm_ppp
R = 0.67231, p-value = 1.194e-11
alternative hypothesis: two-sided